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A detailed numerical study was conducted on the dynamics and thermal response of inert, spherical particles 
in strained, laminar, premixed hydrogen/air flames. The modeling included the solution of the steady 
conservation equations for both the gas and particle phases along and around the stagnation streamline of an 
opposed-jet configuration, and the use of detailed descriptions of chemical kinetics and molecular transport. 
For the gas phase, the equations of mass, momentum, energy, and species are considered, while for the particle 
phase, the model is based on conservation equations of the particle momentum balance in the axial and radial 
direction, the particle number density, and the particle thermal energy equation. The particle momentum 
equation includes the forces as induced by drag, thermophoresis, and gravity. The particle thermal energy 
equation includes the convective/conductive heat exchange between the two phases, as well as radiation 
emission and absorption by the particle. A one-point continuation method is also included in the code thar 
allows for the description of turning points, typical of ignition and extinction behavior. As expected, results 
showed that the particle velocity can be substantially different than the gas phase velocity, especially in the 
presence of large temperature gradients and large strain rates. Large particles were also found to cross the gas 
stagnation plane, stagnate, and eventually reverse as a result of the opposing gas phase velocity. It was also 
shown that the particle number density varies substantially throughout the flowfield, as a result of the straining 
of the flow and the thermal expansion. Finally, for increased values of the particle number density, substantial 
flame cooling to extinction states and modification of the gas phase fluid mechanics were observed. As also 
expected, the effect of gravity was shown to be important for low convective velocities and heavy particles. 
Under such conditions, simulations indicate that the magnitude and direction of the gravitational force can 
substantially affect the profiles of the particle velocity, number density, mass flux, and temperature. © 1999 
by The Combustion Institute 


INTRODUCTION 

Combustion science has been significantly ad- 
vanced during the last 20 years mainly because 
of the evolution and extensive use of advanced 
laser diagnostics along with modeling made 
possible by supercomputers and the associated 
algorithms. Therefore, much has been learned 
in terms of the details of the underlying elemen- 
tary processes of molecular transport and chem- 
ical kinetics as well as their interaction with the 
fluid mechanics. This knowledge, however, has 
been chiefly advanced in the gas phase, in which 
the assumptions of dilute gas can be a very good 
approximation, and many gas properties can be 
accurately determined from first principles. On 
the other hand, there is a wide range of com- 
bustion phenomena of interest in which the gas 
phase interacts in many different ways with a 
liquid or solid phase. While extensive work has 
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been done on sprays, less attention has been 
given to the details of dusty reacting flows. 

Dusty flows are of particular interest for a 
wide range of applications. Particles can be 
present in a gas intentionally or unintentionally, 
and they can be inert or reacting. Inert particles 
in an otherwise reacting gas flow can lead to 
flame cooling and modification of the extinction 
limits of a combustible mixture. Reacting solid 
particles can release substantial amounts of heat 
upon oxidation and can be used either for 
propulsion (e.g., Al, B, and Mg) or power 
generation (coal). Furthermore, accidents can 
occur when a reacting dust accumulates in air 
and which, in the presence of an ignition source, 
can cause explosion. Such explosions can occur 
during lumber milling, in grain elevators, and in 
mine galleries such as the one which occurred at 
the Haswell Colliery and resulted in the famous 
Faraday and Lyell report [1]. 

Particles are also used as seeds in reacting 
flows in order to measure flow velocities by 
using either the laser Doppler velocimetry 
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(LDV) or particle image velocimetry (P1V) 
techniques. The particles that are used in these 
applications must be chemically inert, so that 
they do not alter the chemical composition of 
the reacting stream. Furthermore, they must be 
small and light enough to closely follow the gas 
phase, in order to assure accurate gas-phase 
velocity measurements. 

The addition of solid particles in a flowing gas 
stream can lead to strong couplings between the 
two phases, which can be of dynamic, thermal, 
and chemical nature. The dynamic coupling 
between the two phases is caused by the differ- 
ent inertia possessed by the fluid and solid 
phases, which results in slip between the phases 
that in turn can lead to the development of 
interphasial forces and the modification of the 
velocity fields of both phases. In addition to the 
phase-interaction forces, thermophoretic, cen- 
trifugal, electrostatic, magnetic, and gravita- 
tional field forces can be exerted on the parti- 
cles. Of these, electrostatic and magnetic forces 
can generally be ignored as few combustors 
possess significant electrical and magnetic fields 
and centrifugal forces can be ignored unless 
strong flow rotation exists. On the other hand, 
thermophoretic forces, which are caused by 
steep temperature gradients, can be important 
in reacting dusty flows. The gravitational forces 
are unavoidably present and can play a substan- 
tial role in the overall dynamic response of the 
particle, especially for heavy particles. Indeed, 
existing experiments on dusty flows (e.g., [2]) 
indicate that global flame properties can be 
quite different when they are conducted at 
normal- and micro-gravity. 

To understand the thermal effects of dust on 
combustion, it must be first realized that as the 
thermal capacity of the solid phase is of the 
order of 10 3 times that of the gas phase, it will 
respond more slowly to the temperature 
changes induced by the flame and, because the 
solid and gas phases are thermally coupled, 
slows the thermal response of the gas phase. 
Such temperature modifications can substan- 
tially affect the chemical kinetics of each phase. 
The thermal equilibration process will depend 
on the thermal properties of the two phases. 
Similarly, the emissivity of the solid phase is 
substantially higher compared to that of the gas 


phase so that the radiative transfer from the 
particles can result in significant energy losses. 

The ultimate challenge of this type of re- 
search is to better understand the chemical 
coupling between the two phases. The problem 
is complicated because the actual mechanisms 
leading to particle gasification can differ de- 
pending on the chemical composition of the 
particle. For example, carbon burning initiates 
as a surface reaction that produces CO as a 
byproduct, which then diffuses from the surface 
and undergoes a secondary oxidation, ivlelal 
powders, such as AJ, B, and Mg are first evap- 
orated due to elevated temperatures and their 
vapor undergoes gas phase burning with the 
surrounding [3], The chemical coupling between 
the two phases leads to particle size reduction, 
modification of the gas phase species composi- 
tion, and the elevation of the temperature of 
both phases as a result of the exothermicity of 
both the surface and gas phase reactions. How- 
ever, all particles begin this process as inert 
particles that must first be heated before they 
can themselves participate in the combustion 
process. Thus, one must first understand the 
more basic questions of the interactions of inert 
particles in a combustion environment before 
proceeding to the more complicated chemical 
effects. 

The in-detail understanding of the dynamics 
and structure of heterogeneous (sprays or par- 
ticle) flows can be only advanced by first con- 
sidering simple flow geometries, which (1) can 
be conveniently produced in the laboratory, (2) 
can be simulated with the use of detailed de- 
scription of all the all physico-chemical pro- 
cesses in both phases, and (3) are of relevance 
to practical turbulent reacting flows. Past expe- 
rience in gas phase combustion has shown that 
the opposed-jet, stagnation-type flows are 
among the most meritorious and permit an 
in-depth understanding of the details of the 
pertinent physico-chemical processes. 

A number of computational and experimen- 
tal studies on sprays and particle flows (e.g., 
[4-12]) have been conducted in stagnation-type 
configurations. Numerically, the need for a hy- 
brid Eulerian-Lagrangian approach has been 
identified by Continillo and Sirignano [4], and 
the use of such approach has allowed for the 
prediction of the phenomenon of droplet flow 
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reversal [5. 6], which has been observed exper- 
imentally [7, 8]. Gomez and Rosner [9] have 
conducted a detailed study on the particle re- 
sponse in the opposed-jet configuration, and the 
particle thermophoretic diffusivities were deter- 
mined experimentally. Sung, Law, and co-work- 
ers have conducted numerical studies [10, 11] 
on the effect of strain rate and temperature 
gradients on the dynamics of inert particles, as a 
way of understanding potential errors in exper- 
imental LDV data that might arise from ther- 
mophoretic forces that cause the motion of the 
tracer particles to differ from that of the gas. 
These studies have included the detailed de- 
scription of the particle Stokes drag and ther- 
mophoretic forces. Results have shown that 
depending on the particle size, strain rate, and 
temperature gradient, the particle velocity can 
substantially differ from the gas phase velocity, 
which indeed compromises the fidelity of exist- 
ing non-intrusive flow velocity measurement 
techniques. In these studies, the particle phase 
was not thermally coupled with the gas phase 
and was not allowed to affect the dynamics of 
the gas phase. Furthermore, the effect of gravity 
on the particle dynamics was not considered, 
since the focus of the study were the sub-micron 
sized particles, typically used as tracers, for 
which gravitational forces are negligible com- 
pared to fluid drag. 

The behavior of reacting dusty flows can also 
have a strong dynamic and thermal dependence 
on the particle number density, n p , which rep- 
resents the number of solid particles per unit 
volume; such effects were not addressed in the 
previous numerical studies [10, 11], Nonethe- 
less, the values of n p observed in a combustion 
environment will in most cases be very small. 
This can be easily seen by realizing that a 
stoichiometric mixture of volatile particles and 
oxidizer requires delivering masses of the same 
order of each to the flame, but since the density 
of the solid is roughly 3 orders of magnitude 
larger than that of the surrounding gas, the 
volumetric loading will be of the order of 10~ 3 . 
Consequently, it is unlikely that particles will 
interact with one another and modeling infor- 
mation derived from studies on single particles 
in infinite fluids may be used with confidence. 

In view of the foregoing considerations, we 
have undertaken the task to conduct a long- 


term combined experimental and num 
study on the details of reacting dusty 
which the effects of fluid mechanics 
properties, steep temperature gradients 01 
and particle detailed kinetics, and gravj^ 
systematically addressed. Given the com i 3 
that reacting particles introduce, in 0ui . e 
attempt to understand such complex phe no 
ena, the dynamic, thermal, and gravitation 
effects were independently studied for the ca 
of chemically inert particles; the term “j ne 4 
indicates that the particles are neither a ■ 
sumed nor they alter the gas phase chemistwS 
Thus in this paper, a systematic detailed numejv 
ical study is presented for inert particles embed} 
ded in strained, laminar premixed flames, which 
are stabilized in a stagnation-type flow configu- 
ration. Studies of chemically reacting particles 
will follow. 


NUMERICAL APPROACH 

Governing Conservation Equations 

As mentioned in the Introduction, it is assui ? 
that the particle number density is small enouj 
so that the particles are very unlikely to encotjl 
ter one another which greatly simplifies a 
governing equations. In describing the mechai 
ics, we will consider both the solid and the a 
phase. The gas includes both the original of 
dizer and the gaseous species generated by il 
chemical reactions. The general steady consjg 
vation equations will be presented for the 
of chemically inert particles in axisymmel 
coordinates. Subsequently, these equations 
be reduced to a quasi-one-dimensional fondj 
lation appropriate for the description of stags 
tion-type flows. In these equations, the 
scripts “g” and “p” correspond to the gas || 
particle phase respectively. 

The overall mass continuity equation for 
gas phase is: 


1 a(rp,v R ) 


dr 


d< ' p s u s > _ o 
dx 


The gas phase steady momentum equa 
in the axial and radial directions are, r- 
lively: 
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In the above, equations, u and v represent the 
velocities in the axial (x) and radial (r) direc- 
tions respectively, p. is the mixture dynamic 
viscosity, while and n„F PR are the forces 

in the axial and radial directions, respectively, 
exerted upon a unit volume of gas by the 
particles [13] (the expressions for the forces 
exerted on each particle F PX and F PR in the 


axial and radial direction, respectively, will be 
given later in this section). The gravity term has 
been included only in the x-direction, which is 
parallel to the axis of symmetry of the counter- 
flow. 

The energy equation of the gas phase is given 
by 


r PgiW 


37); 

dx 


+ fpgVgCp 





K qj * K K 

r Pg X ^k^pk^kx ^ r Pg X ^V-pk^kr ^ r r X ^k^k^k £?g.rad 


k= I k=J 

+ «p[2 p + (« P - « g )Fpx + (v p - 

In the above equation, A g is the gas phase 
mixture conductivity, c p is the mixture specific 
heat at constant pressure, c pk is the specific heat 
at constant pressure of species k, Y k is the mass 
fraction of species k, fF k is the molecular weight 
of species k, V ^ and F kr is the diffusion 
velocity of species k in the x and r direction 
respectively, h k is the specific enthalpy of for- 
mation of species k, is the molar rate of 
production/destruction of species k resulting 
rom all gas-phase chemical reactions, and K the 
1013 number of species. The term Qg.rad repre- 
sents the gas radiative loss. The term Q p repre- 
ents the heat exchange (both by convection and 
on uction) between the particles and the gas 
P **• The terms (u p - u g )F PX and (v p - 
e pr represent the work per unit time done 
‘ Cach P ar hcle against the F px and F PR forces, 


k-t 

g)R pr] = o (4) 

respectively. For a reactive dusty flow, however, 
the heat release due to chemical reactions and 
the heat transfer because of steep temperature 
gradients dominate the contributions of the 
work terms so that (u p - u g )F PX and (v p - 
Vg)F PR can be considered as negligible. 

The gas radiation term Q s ra d is given at the 
optically thin limit by 

Q g .rad = 4a g o(r g - Tt), (5) 

where u is the Stefan-Boltzmann constant, T u is 
the ambient temperature, and a g is the total 
Planck’s mean absorption coefficient of the gas. 

Finally, the heat transfer term between one 
particle and the gas, Q p , is given for low (<1) 
particle Reynolds number, Re p , by (e.g., [13]): 

Q P = 4 ird p A g (T g - T p ), 


( 6 ) 
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where d p is the particle diameter. 

The conservation equation for the gas phase 
species k is 




+ rP # « 


dr ) 


+ ~ ( r P& Y k^kx) 


+ £ (rp e Y k V kl ) - rW^ = 0 

( 7 ) 

While the presentation of the gas phase equa- 
tions above is an Eulerian one, the presentation 
of the particle equations in the following will be 
in Lagrangian terms. As first noticed by Conti- 
nillo and Sirignano [4], this can be necessary for 
the simultaneous numerical integration of the 
gas and particle phase set of equations. 

For the momentum balance equations for the 
particles, the formulation of Sung et al. [10, 11] 
was the starting point and appropriate modifi- 
cations were introduced. In the absence of 
centrifugal, electrostatic and magnetic forces, 
the complete momentum balance for a single 
particle in the axial direction, is given in a 
Lagrangian frame of reference by: 

j = ^PX = ^SDX + ^TPX + ^GRX> 

( 8 ) 

where m p is the particle mass and is given by 



ones studied herein [9]. The effect of diffusio- 
phoresis was estimated by using the formulation 
recommended by Gome 2 and Rosner [9], given 
that our calculations were conducted for lean 
H 2 /air flames. As species with very different 
molecular weights are present, the diffusio- 
phoresis will be a function of the equivalence 
ratio, d>. It was found that its maximum value is 
1.4 cm/s for the <f> = 0.57 and 0.4 cm/s for the 
leaner <f> = 0.25 hydrogen/air flames, the two 
cases studied herein. Given that these values 
constitute only a minor correction to the other 
types of velocities, the diffusiophoretic effect 
was not further considered. 

In Eq. (8), F SD represents the Stokes drag 
force which owes its existence to the velocity slip 
between the two phases and is given for low 
(<l)fle p by(e.g.,[10, 11]): 

Fsox = -^ P c Up ~“ g) ~ (10) 

This corresponds to the Type III approximation 
of Hjelmfelt and Mockros [14] with the addition 
of a slip correction factor C, which is needed in 
order to modify the Stokes’ law for the sub-/xm 
seeding particles. The Knudsen- Weber expres- 
sion for C and for all Knudsen numbers is (e.g., 
[ 10 , 11 ]): 

C = 1 + Kn[a + P exp (— y /Kn)], (11) 


^ Ppar (^) 

and p par is the density of the particle. 

Equation (8) does not include the terms 
representing forces due to the pressure gradient 
in the fluid, due to the fluid resistance to 
accelerating sphere, and due to the drag associ- 
ated with unsteady motion. As Sung et al. [10, 
11] have argued, these terms can be neglected 
as they are proportional to the gas density, 
which is much smaller than the particle density. 
Equation (8) does not also include other 
phoretic effects such as diffusiophoresis, elec- 
trophoresis, and photophoresis [9], Any effects 
of electrophoresis was not considered as the 
particles are not charged and no electric field is 
present. The effect of photophoresis was also 
ignored as it has been shown to have negligible 
contribution for conditions analogous to the 


where a, |3, and y are constants derived from 
the fitting of the Knudsen- Weber formula to the 
experimental data and are equal to 1.142, 0.558, 
and 0.999, respectively [15]. The particle Knud- 
sen number is defined as Kn = 2l fp /d p , where /fp 
is the mean free path of the gas molecules given 
by n t = <J>Pgf fp C mg (e.g., [16, 17]) where is a 
constant equal to 0.491, and c mg is the mean 
velocity of the gas molecules. Assuming that the 
molecules follow a Maxwell distribution, c mg is 
given by (e.g., [17]): 


c mg ~ (8& B oLTz7g/' ,rm g) 1/2 - (12) 


In Eq. 12, ^boltz is the Boltzmann constant 
and m g is the reduced molecule mass. For a 
multicomponent mixture m g is given by 
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where 

W k 

m g.k = J7 (14) 

7 v Avogadro 

with /V Avogadro being the Avogadro constant. 

As mentioned earlier, Eq. 10 is applicable 
only for very low Re p numbers. It was one of the 
goals of the present study, however, to investi- 
gate cases in which the particles’ Re p is larger 
compared to the studies of Sung et al. [10, 11]. 
Thus, a correction was introduced to Eq. 10 
(e-g., [18]): 


SDX 


= g - (1 + 0.15Re° 687 ), 

{ (15) 


where Re„ is defined as 

R Cp m ?^ (16) 

/i g 

It should be noted in passing, that the relative 
velocity between the gas and particle phase is 
the appropriate characteristic velocity to use in 
the definition of Re p . Results showed that this 
correction could play an important role on the 
particle dynamic behavior as Re p increases. 

The term F T p X represents the thermo- 
phoretic force on a spherical particle due to gas 
phase temperature gradient VT g . In the near- 
continuum limit, F r p X is given from Brock [19]: 


/A \ VT 

-6-n-pg-pgdpC^ + C t Kn j 

f tpx = / V . 

(1 + 3C m Kn)[ 1 + 2 -^ + 2C x Kn I 

(17) 

where i? g = p g /p g , and C m , C s , and C t constants 
with values 1.14, 1.17, and 2.18, respectively 
[20], which assure that the fitting formula covers 
satisfactorily the entire range of Kn numbers. In 
the work of Talbot, Cheng, Schefer, and Willis 

[20] , this is justified by realizing that Eq. 17 
degenerates within 3% to the collisionless limit 
('+•, as Kn = 2l ff Jd p — » ») given by Waldmann 

[21] . The thermophoretic force F T p X is included 
because of its significance in reacting flows, in 
w hich particles with diameters of the order of 
mi crons are flowing against substantial temper- 
a ture gradients. Then, the thermophoretic force 


has a direction pointing from the high gas 
temperatures towards the low gas temperatures 
and opposes the particle motion towards the 
flame (e.g., [9]). 

The term F GRX represents the gravitational 
force on each particle and is given by: 

F G Rx=- m p8> ( 18 ) 

where g is the prevailing gravitational accelera- 
tion. In the present investigation various values 
of g were used, in order to systematically study 
the effect of gravity. Under normal gravity 
conditions, g = +981.0 cm/s 2 for cases in which 
gravity opposes the particle motion (+g case), 
andg = -981.0 cm/s 2 for cases in which gravity 
favors the particle motion (-g case). Baseline 
studies were also performed at zero gravity g = 

0 C3.SC). 

For a vertically oriented stagnation flow con- 
figuration in which all flame properties are 
radially uniform, the complete momentum bal- 
ance for a single particle in the radial direction 
is given in a Lagrangian frame of reference by 



= F p R = F 


SDR, 


where 


* SDR - r 


(19) 


( 20 ) 


In Eq. 20, Re p number corrections similar to the 
ones introduced in Eq. 15 are not required, 
given that around the centerline the velocities 
v p and v g are very small and their difference 
results in low Re p . 

In order to reduce the system of the gas phase 
equations from the axisymmetric to a quasi-one- 
dimensional formulation similarly to Kee et al. 
[22], the quantities n p , u p , and F PX must be 
considered as functions of x only. In the ideal 
case of a strictly planar boundary layer between 
two opposed jets, this formulation describes the 
flow throughout the spatial domain between the 
two jets. In realistic experiments, however, ob- 
taining planar geometry throughout is not pos- 
sible given that "edge” effects result from the 
unavoidable shear, which develops between the 
counterflowing jets and the ambient stagnant or 
coflowing gas. Furthermore, the experimental 
flames always have a minor curvature around 
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the centerline, which is result of the pressure 
gradients which develop in the stagnation flow 
and which can affect the flatness of nozzle exit 
velocity, » cxit . Under such realistic conditions, 
the quasi-one dimensional formulation is a valid 
approximation at least in the immediate vicinity 
of the centerline of the system, as the symmetry 
requirement imposes that radial gradients must 
be zero asymptotically close to the centerline. 
Furthermore, Eqs. 19 and 20 can be combined 
to 

_ ~3rrp. grfp(Pp ~ v g ) 
pp \dr) C 


problem. Thus, the reduced quasi-one-dimen- 
sional gas phase mass continuity becomes 


2 Pg G g 


df^) 

dr 


(23) 


The reduced quasi-one-dimensional radial mo- 
mentum equation for the gas phase becomes: 


/dG,\ 

d 

' dc; 

( ck ) 


^ dT 


37rMgdp(G p - G g ) _ 


= 0. 


(24) 


and by dividing both sides by r 


m r 



r 



— 3l W*p | y p v ej 


or 

m p G p G p =* (G p - G g ). (21) 

Equation 21 was derived based on the fact that 
the radial velocities of both the gas and the 
particle must vary linearly with r, as a result of 
a Taylor’s expansion asymptotically around the 
centerline (small r). The quantities G g and G p 
have been defined as 


r*t = 

r dr 


and 



( 22 ) 


and are functions of x only, i.e., G g = G g (x) 
and G p = G p {x). 

The analysis conducted by Kee et al. [22] 
shows that if all properties are only a function of 
x then by differentiating the ^'-momentum equa- 
tion with respect to r and the r-momentum 
equation with respect to x, the pressure curva- 
ture J % = (1 lr)(dP/dr) emerges as an eigen- 
value of the problem. In Ref. 22, gravity was not 
included. If the gravity term Pgg is included in 
the x-momentum equation, and given that this 
term is only a function of x, at least asymptoti- 
cally around the centerline, it can be easily 
shown that / g remains an eigenvalue of the 


It should be noted that while Eq. 23 is identical 
to the one derived for the gas phase stagnation 
flow [22], Eq. 24 is modified through the addi- 
tion of the last term in the left-hand side, which 
accounts for the gas-particle forces. 

It is important to emphasize that the presence 
of gravity does not affect the u g velocity profile 
along the centerline. This a result of the re- 
quirement that gravity acts only in the x-direc- 
tion and that radial symmetry must prevail at 
least around the centerline. Thus, for a given set 
of boundary conditions for u exit at the two ends 
of the finite domain, G g is determined through 
the radial momentum Eq. 24 while u g is deter- 
mined from the gas phase mass continuity Eq. 
23. Gravity, however, alters the axial variation 
of the hydrodynamic pressure through the 
dP/dx term, which has to be modified at the 
exits of the two nozzles in order to maintain the 
u ex it values imposed as boundary conditions. 
This physically compensates for the presence of 
gravitational forces. This analysis is only valid 
either for the ideal case of infinitely large noz- 
zles resulting in strictly planar geometry, or at 
least around the centerline for the more rele- 
vant case of finite-size nozzles. 

The quasi-one-dimensional energy equation 
for the gas phase becomes 



k dT K 

+ p g 2 nc pk f / u- s £ + 2 

lt = l k = 1 

+ Qg.rad + n pQp ~ 0 - 


(25) 
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In Eq. 25, Q grad is given by Eq. 5, while the 
particle-gas heat transfer was accounted for 
larger Re p ’ s by using a correction factor C NU : 

Qp = 4W p A g (T g - r p )C NU . (26) 

C NU is a function of Re p that reflects the way 
that the Nusselt number increases with Re p and 
was determined as a function of In (Re p ) by 
fitting a 4th-order polynomial to the experimen- 
tal results reported by Clift, Grace, and Weber 
[23]. This correction is essential when studying 
large particles that typically operate at Re p ’s 
larger than one. 

The quasi-one-dimensional conservation 
equation for the gas phase species k is given by 

/ dY k \ d . . . , 

Pg U g^ J ^k w k — kb (^7; 

The formulation of the thermal energy equa- 
tion for a single particle was based on the 
assumption that the particle temperature is 
uniform and is a function of time only. This is a 
valid assumption for small particles as long as 
the solid conductivity is much larger than that of 
the gas. Thus, the particle thermal energy equa- 
tion is given in a Lagrangian frame of reference by 

n, p c par( j ~ Qp + Q p.rad ~ 2g.rad~p — 0- 

(28) 

The term Q p is given by Eq. 26. The term Q p rad 
is the radiative heat loss emitted from the 
particle surface and is given by 

e p ,ad=/lp e p a(P p -Tt), (29) 

where c par is the particle specific heat, e p is the 
particle emissivity, and/l p = i rd p is the surface 
of the (spherical) particle. 

The term <2 g . rad _ p represents the thermal 
energy, which is emitted by the gas phase and 
eventually is absorbed by each particle. In eval- 
uating this term, each particle at any spatial 
location was allowed to absorb heat that is 
radiated from every point in the gas phase. The 
heat absorption was determined by dividing the 
gas into layers parallel to the flame and then 
computing the radiative flux using a configura- 
tl0n factor for exchange between a circular disk, 
re presenting the gas layer, and a sphere, repre- 
Se nting the particle [24]. Subsequently, £) g . rad _ p 


was determined by summing the absorbed en- 
ergy from all gas layers. 

The mass conservation of the particles is 
described in an Eulerian frame of reference by 


2ppU p + 


d <Pp“p) _ Q 
dx 


where p p = ttptfi p . For constant m p : 


2« p U p + 


d(n £ n E ) 

dx 


(30) 


Eq. 30 can be written in a Lagrangian frame of 
reference as: 


2« p G p + 


dO^lupj) = 

ds 


(31) 


where s is the coordinate which specifies the 
Lagrangian distance traversed by a particle 
since injection; as the particle may reverse di- 
rection during its flight, it can pass through the 
same Eulerian point, x, many times and its 
properties at that point can be multiple-valued. 
However, the total Lagrangian distance s, will 
steadily increase with time. Note that the veloc- 
ity in the direction of s is the absolute value of 
the Eulerian velocity |n p |. The use of s allows for 
the tracking of n p which becomes singular wher- 
ever the particles undergo path reversal. The 
following integration of Eq. 31 results in a 
formula in which the solid properties will be 
single-valued as functions of s: 


n P = n p .mi exp 


2G n + 


d|w P l 

p ds 


dr 


= Vint exp - 2G p dr 


‘p.mj 

• exp 
Given that 


' d|w 


ds 


^ dr 


dirrjJ d/ _ T d|K p | ds 


ds 


ds i/„ 


f l"rl / ||f f \ 
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Description of Singular Points Through One- 
Point Continuation 


then 

^ “ 1 ‘ 2G p d/j, (32) 
' \ Jo / 

where /? p in j and u p in j are the particle injection 
number density and velocity, i.e., their respec- 
tive values at the nozzle exit. 

A analogous approach for the determination 
of rt p for droplets across singular reversal points 
has also been taken successfully by Gutheil and 
Sirignano [6], who studied the droplet behavior 
by dividing the computational domain in 
“sheets,” with each sheet corresponding to a 
droplet path between reversal points. Further- 
more, a rigorous mathematical explanation has 
been provided, which also proves that this sin- 
gularity is an integrable one. 

It is essential that an important clarification is 
made on the definition of the number density 
n p , in cases where there are length scales in the 
flow direction (i.e., the flame thickness), which 
are of the same size as the particles. The 
number density is technically defined as the 
number of particles counted in a volume, in the 
limit where the volume goes to zero. Presum- 
ably, that volume cannot get smaller than a 
particle diameter, but does that mean that n p is 
meaningless over such scales? The answer is 
largely one of definition of the averaging vol- 
ume and is not a problem for the present study 
where the conditions are steady at each Eule- 
rian point. Under such a condition, the extent of 
the averaging in the direction of the flow (the 
direction for which the flame thickness is a 
relevant scale) can be defined as the flow veloc- 
ity multiplied by a short period of time and can 
take that limit as the time goes to zero. Note 
that this is not a redefinition of n p as, because 
the system is statistically stationary, it is exactly 
the same as making a volume average. Further- 
more, it should be clear that when n p is under- 
stood in this sense, there is no change in the 
meaning of any of the mass conservation equa- 
tions presented above, which is the critical issue. 

In summary, the gas phase is governed by 
Eqs. 23, 24, 25, and 27 while the particle phase 
is governed by Eqs. 8, 21, 28, and 32. 


■pjnjl 


‘p.inj 


exp 


The injection of particles in flames can result in 
substantial modification of the properties of the 
gas phase; for example, the large heat capacity 
of the particles dictates that they heat more 
slowly than the surrounding gas. If the injected 
number density n p inj is large enough, this cool- 
ing may lead to flame extinction. As a result, a 
turning point behavior can be expected [251 
under certain conditions in a flame temperature 
vs n p in j diagram. Such turning-point behavior 
cannot be determined through varying n p in j, as 
the extinction point is singular. This can be 
resolved by following the continuation approach 
of Nishioka et al. [26]. Thus, a one-point con- 
tinuation approach was implemented by impos- 
ing a predetermined gas-phase temperature re- 
duction at one point within the flow field, so that 
the n p inj would become part of the solution, and 
its respective injection (initial) value was re- 
moved. The internal point was chosen to be the 
location at which the temperature has maximum 
slope, following the recommendations of Nish- 
ioka et al. [26]. In the present study the one- 
point continuation approach was applied for the 
determination of n p i „j. However, this approach 
can be easily extended to study extinction 
and/or ignition as functions of the particle tem- 
perature, diameter, or any other physical prop- 
erty of the solid phase. 

Physical and Chemical Properties 

In the present investigation, the response of 
inert A1 2 0 3 particles were studied as a model 
case in premixed, atmospheric, lean H 2 /air 
flames by including detailed descriptions of 
chemical kinetics and molecular transport. 

The hydrogen kinetics were taken from the 
GRI 2.1 kinetic mechanism [27], and the mo- 
lecular properties pertinent to the calculation of 
transport coefficients were taken from the San- 
dia Transport Package code [28]. 

For the calculation of the gas phase radiative 
heat transfer, the total Planck’s mean absorp- 
tion coefficient, a g , was calculated based on the 
H 2 0. The Planck mean absorption coefficient 
a HX) for H,0 is given for optically thin condi- 
tions by Tien [29] and Hubbard and Tien [30] as 
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functions of temperature. Then, the total mean 
absorption coefficient, a g is given by 

a g = °h,qPh,o> (33) 

where /» H ,o < s the partial pressure of H 2 0. More 
details on the formulation of the gas radiation 
can be found in [31, 32]. 

The physical properties for the A1 2 0 3 parti- 
cles were taken from Rosenhow, Hartnett, and 
Ganic [33], The particle heat capacity, c par , and 
thermal conductivity, A p , were fitted as functions 
of the particle temperature by using 5th-order 
polynomials. The density of the A1 2 0 3 was 
taken as p par = 3.97 g/cm 3 and the particle 
emissivity as e p = 0.25 [33], 

Method of Solution 

The steady conservation equations for both the 
gas and particle phases were simultaneously 
solved by modifying an existing stagnation-flow 
code [34], which had been used to solve steady 
and unsteady premixed and non-premixed 
counterflowing flames. 

The gas phase equations were solved in a 
manner similar to the Sandia flame codes by 
using a damped Newton method. In terms of 
boundary conditions, the reactant, velocity, 
temperature (always 300 K), and concentrations 
at the nozzles’ exits were specified. Further- 
more, plug flow conditions were assumed at the 
nozzle exits. More specifically, the gas phase 
velocity gradients were set equal to zero, al- 
though relaxation of such conditions is not 
expected to alter the underlying physics of the 
problem. The code was integrated into the 
CHEMKJN [35] and Transport [28] subroutine 
Packages. 

The particle conservation Eqs. 8, 28, and 31 
are first-order differential equations and can be 
marched in time (space) by using only one 
initial (boundary) condition, respectively. The 
radial momentum Eq. 21, is algebraic so that G p 
can be computed directly without any initial or 
boundary conditions. As first noticed by Conti- 
n >Ho and Sirignano [4], the time-marching, La- 
8 ra ngian approach for solving Eqs. 8, 28, and 31 
' S P re ferable given that, under certain condi- 
^° ns ’ tJle space-marching, Eulerian approach 
ot handle singular behaviors in space which 


may not be singular in time. This happens when 
the particles undergo flow reversal. This is pos- 
sible when a particle’s inertia forces it to cross 
the gas stagnation plane (GSP), and penetrate 
into a counterflowing gas stream, which eventu- 
ally forces stagnation and reversal of the parti- 
cle stream. 

For the Lagrangian solution of the particle 
equations, the velocity and temperature of the 
particles were set to be equal to these of the gas 
phase at the nozzle exits (t — 0). Relaxation of 
these assumptions is, again, not expected to 
alter the physics of the problem. The Lagrang- 
ian, time marching was subsequently conducted 
with variable time steps, At. It was found that it 
was beneficial to adjust the time step At so that 
the spatial locations resulting from the Lagrang- 
ian marching would coincide with the spatial 
locations as determined by the Eulerian solu- 
tion procedure for the gas phase. Thus, the need 
for spatial interpolations, which inherently in- 
troduce inaccuracies, was eliminated. 

During the Lagrangian time-stepping, the 
particle axial momentum Eq. 8 was first ad- 
vanced to a new time step so that the new u p was 
determined. During the solution of Eq. 8, an 
iterative procedure was followed which would 
converge to the correct value of « p starting from 
a reasonable initial guess. This was needed 
because of the existence of the Re p number 
term in the Stokes drag force as given by Eq. 15. 
Having determined the new u p value, the parti- 
cle thermal energy Eq. 28 was then advanced for 
the determination of the value of T p at the new 
time step. Similarly to the particle axial momen- 
tum equation, an iteration scheme was also 
needed given that the fourth power of T p is 
involved in the Q p , rad term, and a closed form 
solution cannot be directly obtained. The G p 
was determined directly at all times from the 
algebraic Eq. 21. Subsequently, the particle 
number density Eq. 31 was advanced through 
the analytical expression, Eq. 32. It should be 
noted that while the behavior of Eq. 32 is 
singular at the locations of particle flow reversal 
at which u p = 0, its integration just before and 
just after the singular point can be conveniently 
done as G p is finite throughout the particle's 
Lagrangian path. 
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RESULTS AND DISCUSSION 

The interactions between the inert particles and 
a reacting gas phase can be both dynamic and 
thermal, and they are controlled by a number of 
parameters which have to be varied indepen- 
dently. These parameters are the gas phase 
injection velocity (which controls the strain 
rate) and chemical composition, and the parti- 
cle injection velocity, diameter, and number 
density. Since the present study was conducted 
for A1 2 0 3 particles only, the physical properties 
of the particles were fixed. Otherwise, these 
properties can be also independently varied. 

The numerical simulations were conducted 
for opposed-jet, atmospheric, laminar premixed 
H 2 /air flames, with a nozzle separation distance 
of 1.4 cm (left nozzle located at x = -0.7 cm 
and right nozzle atx = +0.7 cm) and for 4> = 
0.57 and 0.25. The nozzle exit velocities, u exit , of 
the gas and particle phases were identical for 
both nozzles and varied from as low as 14 cm/s 
to as high as 800 cm/s. 

The twin-burner assembly was considered to 
be vertical so that the gravitational forces would 
act along the direction of the system centerline. 
Given the unidirectional nature of gravity, its 
effect on the particle dynamics was considered 
for three cases. The first case was that of ( +g) 
in which the particles were injected from the 
bottom nozzle with the gravity opposing the 
particle axial motion. The second case is that of 
( -g) in which the particles were injected from 
the top nozzle with the gravity favoring the 
particle axial motion. The third case is that of 
zero gravity (0-g). The (+g) and ( -g ) condi- 
tions can be easily produced in the laboratory 
under normal gravity conditions. The (0-g) con- 
ditions can and will be produced in special 
NASA facilities such as a drop tower and/or a 
parabolic-trajectory aircraft. 

In all figures which follow, the particles were 
injected from the ieft nozzle at a spatial location 
x = -0.7 cm with a direction from left to right. 
By injecting particles from one nozzle only, an 
asymmetry was imposed to the system. This was 
intentional as it allows for the possibility of 
particle penetration into the opposing side of 
the GSP which leads to a variety of dynamic and 
thermal phenomena. In the previous studies of 
Sung et al. [10, 11] only small particles were 


studied, which nearly follow the gas phase and 
always stagnate closely to the GSP. It should be 
finally noted, that in the present calculations 
two identical premixed flames are symmetrically 
established, one on each side of the GSP. 

In the present study, the particle diameters 
were varied from 0.3 to 100 pm which is the range 
of interest for our scheduled experiments. Finally, 
the injection number density at the nozzle exit, 
n p inj , was varied from values as low as 10 particles/ 
cm 3 to values high enough to cool the flame and 
eventually cause local and global extinction. 

Dynamic Effects on the Particles 

The particles’ dynamic response is chiefly con- 
trolled by its momentum through Eqs. 8 and 21. 
Equation 8 can be easily transformed to an 
Eulerian frame of reference, and by dividing 
both sides by ujn p and by realizing that m p is 
proportional to (d p ) 3 it can be shown that the 
spatial variation of u p is given by 


du p 

dx 


, v -^2 

u p d\ {u * U * ) + u p dl 


XL , l_ 
r g « p ’ 


(34) 


where A , and A 2 are combinations of various 
parameters contained in Eqs. 15 and 17. A 
careful inspection of Eq. 34 reveals that as d p 
increases the contribution of the thermo- 
phoretic force is reduced. The Stokes drag 
remains important, as a larger d p implies larger 
particle inertia, which can lead to larger velocity 
difference (w g - u p ) and which partially offsets 
the effect of the diameter in the denominator. 
Equation 34 also reveals that as u p increases, 
the thermophoretic contribution is also re- 
duced. Similarly to the previous case, an arbi- 
trarily large u p can result in large slip velocities 
(w g - Up). Thus, the effect of thermophoresis is 
expected to be important for small particles and 
low particle velocities. A comparison between the 
Stokes drag and the gravity term reveals that the 
gravitational force will be important for large 
particles and low particle velocities, as expected. 

Stokes drag and thermophoretic effects 

The effects of Stokes drag and thermophoresis 
are demonstrated in Fig. 1, which depicts the u p 
profiles of 0.3 and 5.0 pm particles respectively 
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Fig. 1. Gas phase and particle axial velocities profiles for a 
<(> = 0.25 H 2 /air flame with u CJlj , = 100 cm/s. /i pJn j = 10 
particles/cm’, and (a) d p = 0.3 pm and (b) d p = 5.0 pm with 
and without the inclusion of thermophoresis. Simulations 
also included the contributions of Stokes drag and normal 
gravity (+£). 

for the ct> = 0.25 flame. Results indicate that for 
d p = 0.3 jxm, the particles follow closely the gas 
phase in the hydrodynamic zone and that inside 
the flame zone a substantial velocity difference 
develops between the two phases. This differ- 
ence was found to be larger for the weaker = 
0-25 flame in which the extent of the thermal 
expansion and the resulting gas phase velocities 
m this expansion zone are lower compared to 
the stronger <f> = 0.57 flame. For the larger 5.0 
Mm particles, it can be seen that the increased 
inertia of the particles results in a substantial 
difference between n p and u g even in the decel- 
erating hydrodynamic zone. Calculations were 
also conducted without accounting for the effect 
°f the thermophoresis and the results are also 
shown in Fig. 1. It is apparent that thermo- 
Phoresis is chiefly responsible for the discrep- 
anc y between w g and n p within the thermal 



Fig. 2. Gas phase and particle axial velocities profiles for a 
<f> = 0.57 H 2 /air flame with n p inj = 10 particles/cm 3 , for (a) 
“exii = 400 cm/s and d p = 10, 20, and 50 pm and (b) u ni , = 
800 cm/s and d p = 10 and 20 pm. Simulations included the 
contributions of Stokes drag, thermophoresis, and normal 
gravity (+g). 

expansion zone and that by neglecting F X p X , « p 
closely follows u g . These observations are con- 
sistent with similar previous studies (e.g., [10, 

H])- 

Typically, both the 0.3 and 5.0 pm particles 
used to generate the data in Fig. 1 possess small 
particle inertia as is desirable for tracer particles 
used in LDV and PIV measurements, and Fig. 1 
shows that these particles can closely follow the 
gas phase and, in particular, that they reach zero 
velocity at the GSP. By further increasing the 
particle inertia, which may be accomplished 
either by increasing the injection velocity and/or 
the particle diameter, it is possible that the 
particles will penetrate the GSP of the gas phase 
and stagnate at a different location, defined as 
the particle stagnation plane (PSP). This is 
illustrated in Fig. 2 for the <5 = 0.57 flame for 
w exii = 400 and 800 cm/s, respectively, and for 
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various particle sizes. Figure 2a depicts that for 
“exit = 4 00 cm/s, the 10 and 20 /xm particles 
penetrate the GSP, and stagnate within the 
region of the opposing jet emerging from the 
upper burner. This leaves them with zero veloc- 
ity in a region of reversed gas flow. Conse- 
quently, the particles reverse direction and 
again cross the GSP where they may stagnate 
again and undergo a subsequent reversal. This 
results in an oscillating motion, but one that is 
strongly damped, so that the particles eventually 
stagnate on the GSP. This oscillatory behavior 
has been previously observed for droplets [5, 6, 
8] and particles [7] in stagnation flow configu- 
rations. The results of Fig. 2a indicate that the 
20 ixm particles penetrate deeper into the op- 
posing jet regime than do the 10 /xm particles 
because of their increased inertia, while the 50 
ix m particles have a high enough inertia to reach 
all the way to the opposing nozzle exit. The 
results of Fig. 2b for the higher u exJt = 800 cm/s 
indicate that the 20 /x m particles reach much 
deeper into the opposing jet regime compared 
to the u exit = 400 cm/s case, and they undergo 
the first reversal at a location very close to the 
exit of the opposing nozzle. Note that “kinks” 
can be observed in the 20 /xm curves for both 
injection velocities; for the u exit = 400 cm/s line, 
the kink can be observed on the first reversal 
point on the right hand side of the GSP, while 
for the u exit = 800 cm/s line, the kink can be 
seen at the second reversal point on the left side 
of the GSP. These kinks correspond to the 
points where the particles pass through the 
thermal expansion zones which surround the 
flames on either side of the GSP. Notice that, at 
those points, the gas undergoes a rapid acceler- 
ation in response to the intense heating. Thus, if 
the particles can pass through the flames, they 
enter into a regime of reduced gas velocity 
within which they decelerate much more slowly. 
The kinks observed in Fig. 2 are a result of this 
rapid change in deceleration rate. 

Number density effects 

The variation of the particle number density, n p , 
was also studied for a wide range of conditions. 
These studies show certain peculiarities espe- 
cially around the PSP where the axial term of 
Eq. 31 becomes singular. Our method of dealing 



Fig. 3. Profiles of scaled by its injection value partii 
number density for a <f> = 0.57 H,/air flame with ti„„ = 1 
cm/s, ti p inj = 10 particles/cm 3 , and d p = 0.3, 5.0, 20. 50. ai 
100 urn. Simulations included the contributions of Stok 
drag, thermophoresis, and normal gravity ( +g). Temper 
ture profile (solid line) indicates flame location. 

with the singular behavior was discussed in th 
Numerical Approach section. 

The spatial variation of n p , scaled by n p inp i 
shown in Fig. 3 for the rf> = 0.57 flame, « exjt : 
114 cm/s, and particle diameters 0.3, 5.0, 20, 5C 
and 100 /xm. Note that only the region to th< 
left of the GSP is shown in this figure. For thi 
0.3-20 ix m particles, it can be seen that n p i; 
characterized by three distinct slopes. In tht 
hydrodynamic region of the flow, it is constan 
or decreases slowly. It then undergoes a rapic 
reduction within the thermal expansion regior 
and subsequently, this reduction becomes 
milder. The rapid reduction in the thermal 
expansion region can be understood as the gas 
stream undergoes a rapid acceleration there, 
which induces an acceleration of the particle 
stream; the reduction is then explained as Eq. 
(32) indicates that np//j pinj varies inversely to 
Up. However, the inertia of the 50 and 100 /xm 
particles is large enough to make them immune 
to ail of these effects, and their velocities seem 
to decrease with a nearly-constant slope, which 
is somewhat affected by the thermal expansion. 
The results of Fig. 3 also indicate that the slope 
of n p in the hydrodynamic zone has a non- 
monotonic dependence on d p which can be 
attributed to the competition between the G p 
and d(u p )/dr terms. More specifically, it be- 
comes steeper as d p increases from 0.3 to 20 
/xm, and it becomes milder as d p further in- 
creases to 50 and 100 /xm. 
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Fig. 4. Profiles of scaled by its injection value particle 
number density for a <f> = 0.57 H 2 /air flame with „ =114 
cm/s, n p j„j = 10 particles/cm 3 , and d p = 50 p.m. Simulations 
included the contributions of Stokes drag, thermophoresis, 
and normal gravity ( +g). Arrows indicate the direction of 
the particle motion. Temperature profile (solid line) indi- 
cates flame location. 


The n p variations can be explained through 
Eq. 32 which indicates that the n p is affected by 
both G p and d(u p )/dc. More specifically, n p 
tends to decrease as G p increases and as the 
slope d(w p )/dr becomes less negative. Physi- 
cally, a greater G p indicates more intense radial 
transport of particles, while a less negative 
d(u p )/dr slope indicates less tendency of the 
particle phase to be “compressed.” For very fine 
particles, e.g., d p < 1.0 gm, the particles follow 
the gas phase very closely, and in regions in 
which the gas phase is incompressible (constant 
density) the particle phase should be also “in- 
compressible,” namely, the n p should be con- 
stant. This is reproduced by our simulations in 
% 3 for the 0.3 pm particles whose n p is 
constant except in the zone, which is character- 
ed by rapid thermal expansion. For larger 
particles that do not follow the gas phase 
closely, the particle phase can have an apparent 
compressibility” (i.e., n p changes) even in re- 
gimes of constant gas phase density. This is 
shown in Fig. 3 for particles equal to or larger 
l han 5 ^m, whose n p is reduced even in the 
hydrodynamic zone of the gas phase. 

The effect of particle reversal on n p is shown 
m Fig. 4 f or t he ^ _ o 57 fj ame with w cxiI = 1 14 
Cm/S ’ an d d p = 50 pm, conditions which result 
several motion reversals for the particles. It is 
PParent that the particle reversal leads to a 


Fig. 5. Gas phase and particle axial velocities profiles for a 
<j> = 0.25 Hi/air flame with u cxl , = 14 cm/s, n pinj = 10 
particles/cm 3 , and d p = 20 i±m for (+g). (0 -g). and (-g) 
conditions. Simulations included the contributions of Stokes 
drag and thermophoresis. 

singular behavior for n p , which is observed to 
increase substantially in the vicinity of the 
PSP’s. Mathematically, this is supported by Eqs. 
30, 31, and 32. At the PSP’s, the spatial deriva- 
tive of Up becomes large and can not be bal- 
anced by the finite values of G p , which is 
independently determined from the radial mo- 
mentum equation of the particle. The results of 
Fig. 4 also indicate that as the particles reverse, 
n p rapidly drops and that continues in each 
subsequent reversal until the n p reaches negli- 
gible values. Note that this singular behavior is 
partially a result of the assumption that particles 
do not interact with one another. Obviously, the 
number density is not truly singular but is 
physically limited by the maximum possible 
particle packing. 

Gravitational effects 

The discussion surrounding Eq. 34 indicates 
that the effect of gravity on the particle dynam- 
ics will be substantial for low convective veloc- 
ities and heavy particles. Simulations were con- 
ducted at (+g), (-g), and ( 0 -g) for the <f> = 
0.25 flame with £< exj , = 14 cm/s and d p = 5, 10, 
20, 50, and 100 /xm. Selected results for the 
particle velocity u p are shown in Figs. 5 and 6 
for d p = 20 and 100 junh respectively. The 
results for small (i.e., ^5 /urn) particles are not 
shown as they generally follow the fluid motion 
and show only a small gravitational effect, only 
noticeable in the hydrodynamic zone. The devi- 
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Fig. 6. Gas phase and particle axial velocities profiles for a 
<b = 0.25 H ; /air flame with = 14 cm/s, /t pi „ t = 10 
particles/cm 3 , and d p = 100 p.m for (+g), (0 -g). and (-g) 
conditions. Simulations included the contributions of Stokes 
drag and thermophoresis. Temperature profile (solid line) 
indicates flame location. 


ations disappears in the reaction zone as the 
increased drag forces, (resulting from the higher 
gas velocities that accompany the thermal ex- 
pansion) overwhelm the gravitational forces. 
The gravitational effect becomes more signifi- 
cant for the 10 and 20 jam particles. For exam- 
ple, Fig. 5 (20 gm particles) depicts that (+g) 
conditions result in a negative d(u p )/dx velocity 
gradient at the nozzle exit; the magnitude of this 
change can be expected to increase with d p as 
the gravitational force increases as (<f p ) 3 while 
the drag forces increase proportional to d p . 
Thus, the (-g) condition results in a positive 
d(u p )/dc gradient at the nozzle exit which again 
increases with d p . The results for the (0-g) 
conditions indicate that u p develops initially 
under the influence of the Stokes drag, as 
expected. 

As also expected, the gravity effect on u p was 
found to be even more dramatic for the larger 
particle diameters of 50 and 100 (Fig. 6) jxm, 
respectively. For both sizes under ( -t-g) condi- 
tions, the particles stagnate at a short distance 
(=1 mm for the present conditions) from the 
nozzle exit, before even reach the GSP, reverse, 
and eventually re-enter the nozzle. For (— g) 
and (0-g) conditions, however, the picture is 
different. Figure 6 shows that, at (0-g) the 100 
/x m particles penetrate deep into the flowfield, 
cross the GSP, stagnate inside the opposing jet 
regime, and subsequently reverse several times 
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Fig. 7. For a i = 0.25 FF/air flame with u exiI = 14 cm/s, 
n p mi = 10 particles/cm 3 , and d p = 20 jtm (a) spatial 
variation of the particle number density scaled by its injec- 
tion value (b) spatial variation of the product of particle 
number density times particle axial velocity scaled by its 
injection value for (+g), (0-g), and (-g) conditions. 
Simulations included the contributions of Stokes drag and 
thermophoresis. Temperature (solid line) profile indicates 
flame location. 


until they reach an equilibrium state. Figure 6 
also depicts that the large 100 /xm particles 
under the favorable (-g) conditions, possess 
high enough inertia to shoot through the GSP 
and reach the opposing nozzle. 

A direct consequence of the gravity effect is 
the potentially large modification in the n p 
distribution, especially for large particles. This 
can be especially seen in Fig. 7a for the 20 pun 
particles, where the substantial discrepancy of 
the n p values between the (+g), (~g), and 
(0-g) conditions can be easily seen. The (+g) 
conditions result in particle accumulation at the 
vicinity of the nozzle exit, while the ( -g) con- 
ditions result in a substantial reduction of n p in 
the same area. This discrepancy was found to be 
even more apparent for the 50 and 100 jx m 
particles. 

These observations can partially explain the 
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» substantia] differences which have been ob- 
» served for the burning rates of reacting dusty 
» flows under normal- and micro-gravity Condi- 
's tions (e.g., [2]). The presence of gravity can 

> jj affect the dynamic response of the particles and 

> | substantially modify the spatial distribution of 

> | the particle mass flux. This can be seen in Fig. 

> ^ 7b, which shows the spatial variation of the 

i product n p *u p , (which is directly proportional 
i to the particle mass flux) scaled by its injection 
i value, for (+g), (—g), and (0-g) conditions. 
» For the case of reacting particles, this is trans- 
» lated to a substantial modification of the rate at 
» ^ which (solid) fuel is supplied to the reaction 
■ s zone and thus one can expect the burning rate 
i S to be significantly affected by the gravitational 
. e orientation. 

,t> 

r" 

Thermal Effects on the Particles 



Fig. 8. Gas phase and particle temperatures profiles for a <t> 
= 0.57 Hi/air flame with u cxiI - 114 cm/s, » pinj = 10 
particles/cm 3 , and d p = 0.3, 20, 50, and 100 jxm. Simulations 
included the contributions of Stokes drag, thermophoresis, 
and normal gravity (+g). 


. The particle thermal response is controlled by 
its thermal energy equation (Eq. 28). When 
m/s, transformed to an Eulerian frame of reference, 
atial it can be shown that the spatial variation of the 
tide P art ’ c * e temperature T p is given by 


“»■ j— ■ - 72 f — 7 , too ) 

and w p^p^parPpar u pdp c parPpar 

ates where A 3 and A 4 are combinations of the 
various parameters contained in Eqs. 26, 28, 
and 29. The convective-conductive terms have 
e g keen absorbed in the A 3 parameter, while the 
•les Station terms have been absorbed into they4 4 
ess parameter. Equation 35 indicates that the spa- 
,$p tial variation of T p is inversely proportional to 
the particle velocity, density, and specific heat, 
t is Thus ’ the higher the values of u p , p par , and c par , 
n the slower is the heating rate of the particle. 
hj s Changing the particle diameter d p has a very 
rm different effect on the convective/conductive 
of and the radiative thermal interactions between 
nd f^e two phases. For the convective/conductive 
g ) ln teraction, the particle heating rate is inversely 

he P r °portional to d p indicating the physical im- 
, n - Portance of the ratio between the characteristic 
in kngth for the heat conduction ( d p ) between the 
be tw ° phases and the particle volume (~d p ), 
itn "'itich is a measure of the particle mass and thus 
'he particle's thermal capacity. For the radiative 
he ""fraction, the heat transfer rate is inversely 


proportional to d p indicating the physical im- 
portance of the ratio between the particle sur- 
face (—d p ), which affects the total radiative 
energy transfer, and the particle’s thermal ca- 
pacity which is again proportional to ( ~d p ). 

The effect of d p on the particle heating is 
shown in Fig. 8 for the <j> = 0.57 flame with n exit 
= 114 cm/s andz/ p = 0.3, 20, 50, and 100'pm. It 
is apparent that the smaller, 0.3 /cm particles 
heat up very quickly and follow the gas phase 
temperature closely. For the larger particles, a 
hysteresis exists, which becomes larger as d p 
increases. For the 50 p m particles which un- 
dergo several reversals, it is of interest to note 
that during these reversals around the GSP 
(a: = 0) the particles eventually are heated up 
to temperatures which are close to those of the 
gas phase. 

The effect of particle reversal on the particle 
heating can be seen more clearly in Fig. 9 for 
the <t> = 0.57 flame, with w exjt = 400 cm/s and d p 
= 20 pm. The particles are heated slowly while 
on the left side of the GSP (.v = 0). After they 
penetrate the GSP they reach a local maximum 
in temperature, and are subsequently rapidly 
cooled to their initial injection temperature just 
as they reach the first PSP. This temperature 
reduction is quite rapid as the particles undergo 
intense convective cooling by the opposing jet 
once they penetrate the GSP. Upon reversal, 
the particles are heated, relatively quickly, back 
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Fig. 9. Gas phase and particle temperatures profiles for a <f> 
= 0.57 H : /air flame with u cxit = 400 cm/s, n pjnj = 10 
particles/cm\ and d p = 20 p.m. Simulations included the 
contributions of Stokes drag, thermophoresis, and normal 
gravity ( +g). Arrows indicate the direction of the particle 
motion. 


to temperatures close to those of the gas phase, 
and as they undergo further reversals around 
the GSP, they eventually reach a maximum 
temperature that is close to the maximum gas 
temperature. Note, that during the first reversal 
path towards the flame following the first PSP 
( x =» +0.22 cm), the particles are heated faster 
as compared to when the particles first enter the 
preheat zone ( x = -0.19 cm). This occurs 
because the particle velocity fields are very 
different in these two regions (as can be seen in 
Fig. 2a). When the particles first enter the 
preheat zone, (0.2 cm < x < 0 cm), their 
velocities are of the order of 150 to 250 cm/s 
while in their first direction reversal (0.12 cm < 
x < 0.22 cm) their velocities range from zero to 
approximately 120 cm/s in magnitude. Conse- 
quently, the returning particles have more time 
in contact with the heated fluid and more 
rapidly rise towards the flame temperature. 

The effect of the gas and particle injection 
velocity, « exit , on the particle heating can be 
further seen in Fig. 10 for the <f> = 0.57 flame 
with d p = 100 /am, and M exit = 114, 400, and 800 
cm/s. Results show that as w exit increases, the 
two flames approach the GSP shrinking the 
region occupied by hot products, but at the 
same time, the maximum flame temperature 
increases, as is typical of these Le < 1 flames 
(e.g., [36]). However, reducing u exit results in 
heating the particles to higher temperatures. 



Spatial Distance, cm 


Fig. 10. Gas phase and particle temperatures profiles for a 
</> = 0.57 H,/air flame with u cxi , = 114, 400, and 800 cm's, 
n p.inj = 10 particles/cm\ and d p = 100 urn. Simulations 
included the contributions of Stokes drag, thermophoresis, 
and normal gravity (+g). 


This indicates that the correspondingly larger 
residence times of the particles in the regions of 
large gas temperature are more important in 
determining the particle temperature than an 
increase of the flame temperature. Thus, for 
very high u exU (e.g., H exit = 800 cm/s), the 
particle temperature is only slightly modified 
from its injection value despite the larger flame 
temperatures. 

The effect of gravity on the particle heating is 
shown in Fig. 11 for large particles and small 
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Fig. 11. Gas phase and particle temperatures profiles for a 
<t> = 0.25 Hj/air flame with u„ xil = 30 cm/s, n pinj = 10 
particles/cm 3 , and d p = 100 jim for ( +g) and (0-g) condi- 
tions. Simulations included the contributions of Stokes drag 
and thermophoresis. Arrows indicate the direction of the 
particle motion. 
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injection velocities (<f> = 0.25. d p = 100 p.m, 
and « exi , = 30 cm/s). For (+g), these heavy 
particles stagnate well before they reach the 
GSP and reverse towards the feed nozzle exit. 
During this excursion, the particles are first 
heated only to about 500 K, as they do not have 
a chance to interact with the hotter regions of 
the flame further downstream. During the re- 
turn path the particles cool more gradually than 
the heating generating a hysteretic behavior 
that reflects the particles’ large thermal capac- 
ity, which delays both their initial heating and 
the subsequent cooling. At (0-g), the particles 
penetrate the GSP and undergo several rever- 
sals before settling on the GSP. During these 
reversals, the particles follow a complex temper- 
ature variation as they sequentially interact with 
hotter and colder gas layers. These results un- 
ambiguously show that gravity can have a po- 
tentially strong effect on the particle thermal 
response as well as the particle dynamics. 

The relative importance of the various terms 
of Eq. 28 on the particle thermal state was also 
analyzed in detail. It was found that the convec- 
tive/conductive term, Q p , strongly dominates 
within the flame zone for all cases studied. 
Away from the flame zone, the Q p term was 
found to be comparable in magnitude to 
Qgrad-p as radiation provides a mechanism for 
long-range interactions with the high tempera- 
ture regions of the flow. The particle radiation 
term, £? p . rad , was found to be maximum at the 
vicinity of maximum T p , as it is expected, but in 
such regions its contribution is substantially 
lower than that of Q p . 

Dynamic and Thermal Effects on the Gas 
Phase 

For low particle number density n p , the dynam- 
ics and thermal response of the gas phase are 
nearly unaffected by the presence of particles. 
This was the case for all the results which were 
shown in Figs. 1 through 11. For these simula- 
tions, the low value of n p jn j = 10 particles/cm 3 
w as chosen in order to isolate the particle 
effects from the gas-phase effects. However, as 
increases, the gas phase can be affected by 
the presence of the particles both dynamically 
and thermally, as should be apparent from Eqs. 
24 and 25, respectively. 


The dynamic coupling results from the forces 
which develop between the two phases and thus, 
can affect the momentum balance of the gas 
phase and the gas phase velocity field. This 
possibility was assessed for the <f> = 0.57 flame, 
with u exjt = 400 cm/s, d p = 50 p, m, and n p inj = 
5800 particles/cm 3 . The simulations were con- 
ducted with and without the particle-gas force 
interaction in Eq. 28. It was found that if the 
particles penetrate the GSP, they can affect u g 
especially in a region just to the right of the 
GSP. Thus, by allowing the particle-gas interac- 
tion, u g becomes lower in magnitude. Physically, 
this is a result of the force exerted between the 
two phases, which reaches a maximum just to 
the right of the GSP where the relative velocity 
between the phases, and thus the interphasial 
drag forces, are at a maximum. As a conse- 
quence, the gas experiences a force of equal 
magnitude and opposite sign to that experi- 
enced by the particles, resulting in a reduction 
of the magnitude of u g . 

The thermal effects on the gas phase was 
found to be significant as n p increases. More 
specifically, by injecting particles at ambient 
temperature, the flames are cooled, and, as 
”p.inj * s gradually increased, near-extinction con- 
ditions were observed. The cooling effect was 
found to be stronger for the left flame of the 
twin-flame assembly, as the left flame interacts 
first with the incoming cold particles (which are 
only injected from the left jet). It was also 
shown, that this cooling effect was different for 
particles with small and large inertia, as the 
particles with large inertia penetrate deeper 
into the flowfield and can interact thermally 
with both jets. 

The temperatures of the gas, T g . and parti- 
cles, 7 p , are shown in Fig. 12 for the cf> = 0.57 
flame, with w exjI = 114 cm/s, d p = 20 /rm. and 
"p.inj = 10 and 24,500 particles/cm 3 . It can be 
seen that for the higher « p jn j, T g is substantially 
reduced for the left flame, while the tempera- 
ture of the right flame is nearly unaffected. For 
the higher n p in j, T p is also reduced as a result of 
the reduction in flame temperature. 

In Figs. 13 and 14 results are shown for the <f> 
= 0.57 flame, with « exjt = 400 cm/s and d p = 50 
(im. For this case, the particle inertia is large 
enough that the particles can reach all the way 
to the opposing nozzle exit and thus, thermally 
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Fig. 12. Gas phase and particle temperatures profiles for a 
<b = 0.57 H ; /air flame with u oxjl = 1 14 cm/s, n pinj = 10 and 
24.500 particles/cm 3 , and </ p = 20 jim. Simulations included 
the contributions of Stokes drag, thermophoresis, and nor- 
mal gravity ( +g). 

affect the entire flowfield. The dependence of 
the maximum flame temperature, T f max , on 
« p . inj can be seen in Fig. 13; for all n pJnj ’s, T f max 
is the maximum temperature of the right flame 
as the left flame undergoes more intense cool- 
ing. As expected, the results of Fig. 13 indicate 
that T f max is first monotonically reduced with 
/ip inj. It was found that the left flame is extin- 
guished for n p in j = 13,450 particles/cm 3 , while a 
“turning-point” behavior is observed at n pin j = 
26,627 particles/cm 3 , which indicates extinction 
of the right flame as well and as a consequence 



Fig. 13. Maximum flame temperature variation with parti- 
cle number density, for a <f> = 0.57 H ,/air flame with u e , u = 
400 cm/s, and d ? = 50 nm. Simulations included the 
contributions of Stokes drag, thermophoresis, and normal 
gravity (+g). 
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Fig. 14. Gas phase and particle temperature profiles for ; 
= 0.57 Fi ,/air flame with d p = 50 /am. n cxjl = 400 cm/s. a 
n p.i n j = 10, 7,010, 13,750, and 25.258 particles, ■'em 3 . Simu 
tions included the contributions of Stokes drag, them 
phoresis, and normal gravity (+g). 


global extinction. This “turning-point" extir 
tion behavior is typical of flames undergoi 
cooling and it could be captured by the or 
point continuation technique that was imp 
mented in this model. 

Figure 14 shows the thermal structures of t 
gas and particle phases for three of the cases 
Fig. 13: forn pin j = 10, 7,010, 13,760, and 25,2 
particles/cm 3 . For n p inj = 10 particles/cm 3 ca 
T g is unaffected by the particles, while for t 
n p jnj = 7,010 particles/cm 3 case, T g is substt 
tially reduced throughout the flame asseml 
with, as would be expected, a more profou 
effect on the left side of the flame assembly tl 
first encounters the particles. T p is also smal 
than the /t pin j = 10 particles/cm 3 case, ag; 
reflecting the reduction in flame temperatu 
Thus, these particles remove thermal enei 
from the high-temperature regions of the j 
phase and transfer that energy to other parts 
the flowfield, preheating the gas flow com 
from the right nozzle. This can be seen in f 
14, as the gas phase temperature of the 
emerging from the right nozzle, starts increas 
well ahead of the flame reflecting the prehe 
ing induced by the heated particles. The n p in 
13,760 particles/cm 3 case describes a state j 
after the left flame has been extinguished, i 
results for both the gas and particle pha 
reflect the existence of only a single flame. 1 
n p jn j = 25,258 particles/cm 3 case describe! 
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state just before the right flame extinction, (i.e., 
when both flames are extinguished). The sub- 
stantial reduction of the temperatures of both 
phases is apparent. 

Finally, it should be noted that for all studies 
on the effect of /i p inj on the gas phase, special 
care was taken to account for the local n p 
modification when particle reversal(s) wa- 
s(were) observed. Thus, the effective n p which 
was used in the gas energy equation, was ob- 
tained by summing over the contributions to n p 
resulting by all reversal paths. 

| CONCLUSIONS 

In the present study, the dynamic and thermal 
interactions between inert A1 2 0 3 panicles and 
strained, laminar, premixed H 2 /air flames were 
numerically investigated in a stagnation flow 
configurauow. In this analysis, the quasi-one- 
dimensional conservation equations that de- 
scribe the ideal stagnation flow for the gas phase 
were solved. These equations include the parti- 
cle-gas phase force interaction in the momen- 
tum equation, and detailed description of chem- 
ical kinetics and molecular transport. A set of 
conservation equations was also developed for 
the particles in ideal stagnation flows. The par- 
ticle momentum balance included the effects of 
thermophoresis and gravity in addition to the 
Stokes drag force. The particle thermal energy 
balance included convective/conductive and ra- 
diative exchanges between the particle and the 
surrounding gas. A one-point continuation ap- 
proach was also implemented, allowing for the 
determination of singular, turning points, de- 
scribing flame extinction. The effects of the gas 
phase stoichiometry, strain rate, particle inertia, 
thermophoresis, and gravity were assessed. 

As expected, results show that while small, 
submicron particles follow the gas phase closely, 
larger particles fail to do so. For the small 
Particles, the effect of thermophoresis was 
found to be significant in regions of large tem- 
perature gradients in the gas phase, an effect 
that was not as important for large particles. 

This is in agreement with results of previous 
studies. Particles possessing large inertia were 
found to penetrate the GSP, stagnate in the 
interflow the opposing jet, and subsequently 
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reverse direction several times around the GSP, 
in agreement with previously observed flow 
reversal for droplets in stagnation-type flows. 
Finally, the results demonstrate that the param- 
eters affecting the particle velocity can also 
directly affect the particle number density which 
varies significantly throughout the flowfield. 

Simulations at (+g), (0-g), and (— g) showed 
that the magnitude and direction of the gravi- 
tational force can have a strong effect on heavy 
particles moving at low velocities. In addition to 
significantly modifying the particle velocity field, 
gravity was also found to be responsible for the 
modification of the particle temperature, num- 
ber density, and mass flux distribution, which 
may partially explain experimentally observed 
effects of gravity on flame propagation in dusty 
reacting flows. 

The thermal effects between the two phases 
were also found to be strong. Results show that 
the particle heating rate is lower for high-inertia 
particles. Under conditions that result in parti- 
cle flow reversal, the particle temperature was 
shown to vary in a highly non-monotonic man- 
ner, by interacting with gas layers of vastly 
different temperatures. 

Finally, simulations with high particle number 
densities showed that the presence of the parti- 
cles can affect both the velocity and tempera- 
ture fields of the gas phase. The velocity field 
was found to be affected because of the forces 
that develop between the two phases. The effect 
on the gas phase temperature is mainly a result 
of the heat transfer between the hot flame and 
the cold particles. As the number density in- 
creases, the flames can be cooled to extinction 
states. Under conditions of high particle inertia, 
the heated particles can penetrate deep into the 
opposing jet, and for very high values of the 
number density can cause an “early” preheating 
of the incoming gas. 
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